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Abstract 

In this paper, we analyze a retrial queueing system with Batch Markovian Arrival 
Processes and two types of customers. The rate of individual repeated attempts from 
the orbit is modulated according to a Markov Modulated Poisson Process. Using the 
theory of multi-dimensional asymptotically quasi-Toeplitz Markov chain, we obtain the 
stability condition and the algorithm for calculating the stationary state distribution 
of the system. Main performance measures are presented. Furthermore, we investigate 
some optimization problems. The algorithm for determining the optimal number of 
guard servers and total servers is elaborated. Finally, this queueing system is applied to 
the cellular wireless network. Numerical results to illustrate the optimization problems 
and the impact of retrial on performance measures are provided. We find that the 
performance measures are mainly affected by the two types of customers’ arrivals and 
service patterns, but the retrial rate plays a less crucial role. 
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1. Introduction 


Most queueing systems assume that the arrival process is a stationary Poisson process. But 
such a process does not characterize the typical features of traffic in modern telecommunication 
networks such as correlation and burstiness. The batch Markovian arrival process (BMAP) is a 
useful and appropriate model for describing such features. The BMAP is a generalization of many 
well-known processes including the Markovian arrival process (MAP), the Markov-modulated Pois¬ 
son process (MMPP), the PH-renewal process and the Poisson process. Furthermore, the BMAP 
preserves the tractable Markovian structure. The or igins o f the BMAP can be traced to the devel¬ 


opment of the versatile Markovian point process by 
were introduced in 


Neut; 


(19791). The currently used denotations 


Neuts ( 19891 ). It is widely employed in research for queueing systems, inventory 


systems, reliability engineering, manufact u ring s ystems, computer communication systems, insur¬ 


ance problems and so on. 


Chakravarthvl (1199911 gave a detai l ed sur vey of the results about the 


queues which the input flow is BMAP. 


Hev man &: Lucan toni (20031 ) presented an application of 


the BMAP for modeling the information flows in the modern telecommunication networks. 

Retrial queues are very good mathematical models for cellular wireless networks, telephone 
switching systems, local area networks under the protocols of random multiple access, computer 
systems for competing to gain service from a central processor unit, etc. Over recent years it has 


be en a r a pid g r owth in t 


by 


Fali n (1999), 


re literature on ret rial q ueues. T 


A rtale io 


(tl99ft 


re reader can refer to the survey papers 


2010b lYang &; Templetonl (|1987u for a review of main results and 


methods. Most of the papers and books are focused on single server retrial queues. Multi-server re¬ 
trial queues are more complicated for research and even the model M/M/c is studied until now only 
algorithmically via different truncation schemes. The multi-server retrial queues with the BMAP 
are still not well investigated in literature. In 


Breuer et al. 


(2002|), the classical BMAP/PH/N re¬ 


trial queue was 
Markov chain. 


analyzec 


He et al 


ry means of discrete-time multi-dimensional asymptotically quasi-Toeplitz 


2000 ) gave th e sta ble co nditio n for the BMAP/PH/S/S+K retrial queue 


with PH-retrial times. Recently, 


Klimenok et 


with impatient customers. Later, 


Kim et al. 


al, 


( 2007 1) analyzed the BMAP/PH/N retrial queue 


Markovian flow of breakdowns. Subsequently, 


(2008) j;tuc 


Kim et al 


ied_ the BMAP/PH/N retrial queue with 
( 201 0) considered_the BMAP/PH/N re¬ 


trial queueing system operating in Markovian random environment. 


Dudin &; Klimenok (2012) 


investigated the BMAP/PH/N retrial queue with Markov modulated retrials. Particularly, for the 
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retrial queues with two types of customers, the single server retrial queues have receivec 
attention in the literature. For example, using the supplementary variable method, 


considerable 


Choi Sz Park 


(_199Q|) dealt with an M\, M 2 /G/I retrial queue with two types of customers and non-preemptive 
resume priority. After that, Choi and Pa rk’s model was extend ed to an M \, M 2 /G \, G 2 /I retrial 


queue with priority customers by 


Falin et al 


(1991a) • Later, 


Choi et al 


5) i 


(|1995i l investigated an 


M/G/l retrial queue with tw o types o f cus tomers and finite capacity. Furthermore, using the 
matrix analytic method, 


Martin & Artaleio 


1993) stu 

died an Mi, Mo/Gi, Gb/l/l retrial queue 

. In addition. 

Choi & Chang 

(1999 

) presented a survey of single server 


retrial queue with two types of customers and priority. Recently, 


Wu & Lian (J2013I) introduced 


negative customers into the retrial queues with two types of customers. However, to the best of our 
knowledge, little work has been done on multi-server retrial queues with two types of customers, 
and no work on BMAP/PH/N retrial queue with two types of customers is found in the queueing 
literature. 

Our model can be successfully studied by means of asymptotically quasi-Toeplitz Markov chains. 
Here, we briefly introduce the continuous-time multi-dimensional asymptotically quasi-Toeplitz 
Markov chains. For a comprehensive and excellent overview of multi-dimensional asymptotically 


quasi-Toeplitz Markov ch ains with discrete and continuous time, rea ders may refer to th e pape r 
ini 


Klimenok &; Dudinl (2006j). The next definition and propositions follow iKlimenok &; Dudinl (2006). 
Let = {h; a f} 5 t > 0, be a regular irreducible continuous time Markov chain with the state space 
Q = {(z,a),a £ Li,i = 0,1, - - - , z°;(z, a),a £ L,i > z 0 }. Enumerate the states of the chain 
t > 0} as follows: the states (z, a) are numerated in ascending order of the component i and for 
fixed z, the states (z, a) are arranged in lexicographic order. We denote the block matrix Q = (Qij) 
as the generator of the chain, where the block Qij is formed by intensities <?(i, a ) ; ( y,b) °f the chain, 
transition from the state (z,a) to the state (j. b). 

Definition 1.1. A regular irreducible continuous time Markov chain > 0} is called asymp¬ 
totically quasi-Toeplitz Markov chain if 


(i) Qi,j = 0 for j < i - 1, z > 0. 

(ii) There exit matrices Y n ,n > 0, such that 

\ n — hill A ? - Qi,i+n—1; ^ — 0, 2, 3, • • • , 

2—>• OO 
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Y\ = lim A t 1 Qi i + I, 

i—> oo 

and the matrix is stochastic. 

(in) The jump chain {(, ni n > 1} of the process {£t,t > 0} is non-periodic, 

where the diagonal matrix Aj defines the diagonal entries of the block —Qu, i > 0. 

Let y(z) = be the matrix generating function characterizing the limiting Markov 

chain. The following two propositions present the ergodicity conditions for the Markov chain 
{£t,t > 0} in terms of the generating function Y(z) based on whether the matrix y(l) is an 
irreducible matrix or not. 

Proposition 1.1. IfY(l) is an irreducible matrix and assume that 

(i) the series J2^=i n Qi,i+n-i^ converges for i = 0,i°, 

(ii) the series Y^=i n Qi,i+n-i^ converges for i > i°, and there exists a positive integer i* > i° 
such that this series converges uniformly in the region i > i*, 

(Hi) the sequence Af 1 , i > 0, has an upper bound, 

then the sufficient condition for the Markov chain [ftf > 0} ergodicity is the fulfillment of the 
inequality 

[det(z/ - Y(z))}' \ z= i > 0. 


If y(l) is a reducible matrix, then it can be presented in a normal form which has the following 
structure: 


f yW(z) 


Y(z) = 


Y^(z) 


\ y^ + 1 Y( z ) y(*+b 2 )(z) 




y(b(z) 

y (z) y( 1+1 \z) / 


where Y^\z), i = 1,1, are irreducible square matrices. 
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Proposition 1.2. IfY( 1) is a reducible matrix of the above form and the conditions of (i)-(iii) of 
Proposition 1.1 are satisfied, then the sufficient condition for the Markov chain {ft A > 0} ergodicity 
is the fulfillment of the inequality 


det (zI-Y^(z)) 


2=1 P 0, 


i = 1,1. 


The rest of the paper is organized as follows. The model description is given in Section 2. 
The Markov chain associated with the system is analyzed in Section [3] where we formulate the 
model as a multi-dimensional asymptotically quasi-Toeplitz Markov chain and obtain the ergodicity 
condition of the Markov chain. The algorithms for calculating the stationary state probabilities 
are elaborated in Section 4. A set of performance measures is demonstrated in Section 5. Two 
optimization problems and an algorithm for calculating the optimal values are provided in Section 

6 . An application for the model under discussion and some numerical examples are shown in Section 

7. Finally, summary of the results is presented in Conclusion. 


2. Model Description 


We consider a multi-server retrial queueing system. Two types of customers (primary customers 
and priority customers) arrive to the system according to two independent BMAPs (denoted by 


BMA 


Neuts 


Pi and BMAPf). The notion of the BMAP and its detailed description were given by 
(1989). We denote the directing process of the BMAP\ by { 17 , f > 0} with the state 


space {1, - - - ,W} and the directing process of the BMAP 2 by {7 t,t > 0} with the state space 
{1, • • • ,V}. The two directing processes behave as two independent irreducible continuous-time 
Markov chains. Denote the matrix sequences (D n : n £ No) and ( E n : n £ No) as the sequences of 
characterizing matrices for the BAIAP\ and BMAP 2 , respectively. Introduce the matrix generating 
functions D(z ) = an d E(z) = y \ z\ < 1. We assume the matrix functions 


D(z) and E(z) satisfy all assumptions of Lucantonil (1991). Therefore, the matrix D( 1) is the 


generator of the process {z g,t > 0} and the matrix E{ 1) is the generator of the process { 7 t,t > 0}. 
Suppose that the two arrival processes BMAPi and BMAP 2 start with the initial phase probability 
distributions aq and a 2 , respectively. Let 6\ and (9 2 be the stationary probability vectors of D( 1) 
and E{ 1), respectively. Then Ai = 6\D'{l)e is the stationary arrival rate of the primary customers 
and A 2 = 0 2 P / (l)e is the stationary arrival rate of the priority customers. Moreover, we define 
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A&i = 9i(—D 0 )e and Xb 2 = 02 (—E o )e as the intensities of group primary arrivals and group priority 
arrivals, respectively. Here e is a column vector of appropriate size with all elements equal to 1. 


The service facility consists of c (c > 2) identical servers which are identical and operate 
independently of each other. All customers have the same independent, identical PH distributed 
service times S which is governed by the irreducible continuous time Markov chain {mt, t > 0} with 
the state space {1, • • • , M}. The transitions of the Markov chain {mt, t > 0}, which do not lead to 
service completion, are defined by the irreducible matrix S of size M x M. The transitions of the 
Markov chain { mt, t > 0}, which lead to service completion, are defined by the vector Sq = — Se. At 
the service beginning epoch, the state of the process {mt, t > 0}, is determined by the probabilistic 
row vector ? of size 1 x M. The service time distribution function has the form B(x) = 1 — ge Sx e. 
The mean rate of service is /j, = [— <^5 , ~ 1 el~ 1 . A more detailed information abo ut th e PH typ e 


distribution can be seen in the books 


\euts (19811). Latouche A Ra m aswami (1999) and Hr (2014). 


When a batch of priority customers arrives at the system and there are several servers being 
idle, the priority customers occupy the corresponding number of the servers. If the number of the 
idle servers is insufficient (or all the servers are busy) the rest of the batch (or all the batch) enters 
the retrial group (called orbit). When a batch of primary customers arrives at the system and meets 
no more than g — 1 (1 < g < c — 1 ) busy servers, the primary customers occupy the corresponding 
number of the servers (If the number of the idle servers is insufficient, the rest of the batch goes 
to the orbit); otherwise, the whole batch of primary customers is blocked and goes to the orbit. 
These customers in orbit are said to be repeated customers and the orbit capacity is assumed 
to be unlimited. For all customers in the orbit, the individual repeated attempts are governed 
by a common Markov Modulated Poisson Process (MMPP). The directing process of the Markov 
Modulated Poisson Process is denoted by {rt, t > 0} which is an irreducible continuous time Markov 
chain with the state space {1, • • • , i?}. Transition intensities of the Markov chain {rt, t > 0}, which 
are not related with customers retrials, are defined by the matrix Tq, and transition intensities, 
which are accompanied by retrial, are described by the diagonal matrix T\ = diagjcri,--- ,o\r}. 
The matrices T = Tq + ^i is the infinitesimal generators of the processes {ry,t > 0}. That is, when 
the process rt stays in the state r, r = 1, • • • ,R, each customer from the orbit makes attempts 
to seek the service at exponentially distributed time intervals with mean l/a r . Hence, the total 
rate of retrials is equal to na r ,a r > 0 , when the process rt is in the state r and the orbit size (the 
number of calls on the orbit) is equal to n, n > 0. Let 9q be the stationary probability vector of 
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T. Then a = # 0 ^ 1 e is the stationary retrial rate of the repeated calls in the orbit. On retrial, an 
orbiting customer obtains service immediately from one of the idle servers if the number of busy 
servers is no more than g — 1; otherwise it will return to the orbit for later retrial. We understand 
that there are c — g servers (called guard servers) opening only for the priority customers. 


For the use in the sequel, let us introduce the following notations: 

* e is a column vector of ones of suitable size. When needed we will identify the dimension of this 
vector with a subindex; 

* I is an identity matrix of appropriate size. When needed we will identify the dimension of this 
matrix with a subindex; 

* Ol is a zero matrix of size L x L; 


* diag{a r ,r = 1 ,L} is a diagonal matrix with diagonal entries a r , r = 1 ,L; 

* (g> and © are the symbols of Kronecker product and sum of matrices, see, e.g. 

i > 1, n®° = 1,n®° = o ; 


Graham (1981); 


1 -1 


-k Q® 1 = In m < 8 > ^ ® I n i-m- 1 , l > 1 , for the matrix Q having n rows. 

m =0 


Our work objectives are to derive the stability condition, to calculate the stationary state 
distribution, to derive the main performance measures of the system and optimize the values of g 
and c. 


3. Analysis 

In this section, we first adopt a layered approach for the description of the state space to 
construct the infinitesimal generator Q of the Markov chain which can represent the model. We 
then derive the ergodicity condition of the Markov chain. 

Let 

* ot be the number of repeated customers in the orbit, at the epoch t, t > 0 , ot > 0 ; 

* bt be the number of busy servers, at the epoch t, t > 0, bt = 0 , c; 

* rt be the state of the directing process of the MMPP of repeated customers, at the epoch t, t > 0, 
n = l ,R; 

* vt be the state of the directing process of the BMAP of primary customers, at the epoch t, t > 0, 

v t = W; 
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* 7 t be the state of the directing process of the BMAP of priority customers, at the epoch t, t > 0, 

It = hV; 

* m)- be the state of the directing process of the service on the j th busy server, at the epoch t, 
t > 0, rrif * = 1, M , j = 1, bt . Here, we assume that the busy servers are numerated in order of their 
occupying, i.e., the server, which begins the service, is appointed the maximal number among all 
busy servers; when some server finishes the service and is released, the servers are correspondingly 
enumerated. 

Then, the process of the model can be described by the regular irreducible continuous time 
multi-dimensional Markov chain: 


f = {£t,t > 0} = | o u b u r t ,v t ,lt-, 


m 


(i) 


, m. 


(M 


t > 0 


with the state space 


o > 0, b = 0, c, r = 1, R, v = 1, W, 7 = 1, V, m^ = 1, M, i = 1, 6 j , 

and the dimension of the state space of the Markov chain , t > 0} is equal to K = RWV^y^w~- 
We note that the dimension K can be very large. E.g., if the array (c, R,W,V, M) is equal to 
(10,2,2,2,2), then K = 16376. 


We suppose that the states of the Markov chain {£t,t > 0} are enumerated in lexicographic 
order. By means of considering the probabilities of the Markov chain transitions during an in¬ 
finitesimal time interval, the infinitesimal generator matrix Q of the Markov chain {£f,t > 0} has 
the following structure: 


Q = 


Qoo 

Q01 

Qo2 

Qo 3 

Q10 

Q11 

Ql 2 

Ql 3 


Q21 

Q 22 

Q23 



Q 32 

Q 33 


... \ 


(3.1) 


V 



where the blocks Q,j of size K x I\ are defined by: 

O, V < 

Irwv^S® 1 , V = 

To © Dq © Eq © S — iT\ © I\yvm l > l' = 

(Qi,i)i,i' = { ToSDoffiToffiS®', /' = 

I R ®(D r ®E r )®I M i®g® r , V = 

Irw © -E7 © I M i © 5® r , V = 

Irw © T r © I M i © ?® r , 


l = 

2 ,c, 

l = 

1 7c, 

0,5 

-1, 

57 c, 


r = 

1,5 “ l . 

r = 

5 + 1 - 


z > 0, 


l' = l + r, r = 1 , c — l, l = g, c — 1 . 


Qi,t-1 - ^ g_l 


Ti <8> /vrv ® s 


Ti (gi Iwvm ® ? 


© ® Iwvm s- 1 ® 


Qiji+fc 


9-1 


C V 

z > 0, k > 1. 


Iff <S> Rfe + t, ® /v ® 


fl ® ^fe+s-1 ® ly m < 


> ? ®(s-i) 


h? ® ^fc+1 ® Ems- 1 ® ^ 
Ir® D k <S> Ivmh 


Irw ® -Efe+c ® <© 


iflvr ® -fc'fe+c-i ® Im ® <> 


(c-l) 


hz ® Rfc ® Ivm c_1 Irw ® -^fc+i ® Im 0-1 ® 

Ir® D k ® Wm c + hw ® Ek® Im c ) 


Remark 3.1. When E n = 0, n = 0,1, ■ ■ ■ and c = g, the present system reduces to a re¬ 
trial BMAP/PH/c queueing system with Markov modulated retrials. It may be noted that the 
results above after putting E n = 0, n = 0,1, • • • and c = g agree with the result presented in 


Dudin & Klimenok 201 
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In the following sections, we analyze the queueing system by using the theory of multi-dimensional 


asymptotically quasi-Toeplitz Markov chains. 

The diagonal blocks of the generator Q are given by 

f To © Do © Eo © S — iT © I\vvm 1 : 1 = 0, — 1, 

\Qii)l,l \ , i>0. 

{ T 0 ®D 0 ®E 0 ®S® 1 , l=g 7 c, 

Denote the diagonal entries of the matrices Do, Eo and S as {—A^,u; = 1 ,W}, {—X^\v = 
1, V} and {— s m ,m = 1 ,M} respectively. Then the diagonal matrix A* in Definition 1.1 is given 
by: 

Ai = A + ifj, (3.2) 

where 

A =diag{diag{a r , r = 1, R} © diag{\^\w = 1, W} 

© diag{\^\ v = 1, V} © [diag{s m , m = 1, M}]® 1 , l = 0, c}, 


Xi — diag{Ti <S> Iwvm 1 ^ ~ 0> c }) 


i= ( Irwv e?Ao Mn 

\ Or\VV YI n=g M n 

Further, introduce the matrices I, T and T^, k > 0, of size K x K: 

O r 


I=1-1= 


'RWV Xln=0 M 


I 


RWV En=9 M n 


T = 


0 

1 

9-1 

9 


Orwv Irwv © ? 

Orwvm Irwvm © g 


' • Irwvmq - 1 © ? 

O RWV Ms 


OrwVM c ) 
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o 


(To )l,V 


V < l - 1,1 = 2,c, 

Irwv <8> S® 1 , l' = l — 1, Z = 1, c, 

T 0 ® £>00^0©^, Z' = Z,Z = 07c, 

Ir <8> (Dr ©Dr) ®/ M i ©xT 7- , Z' = Z + r, r = 1,5 - l,l = 0, g - 1, 

/rw 0 E r ®I M i V = Z + r, r = 3 + 1 — l,c - l, l = 0,g - 1, 

Iru/ <g> E r (g> I M i <g> <^ 8 ’ r , l' = l + r, r = 1, c — Z, Z = 5 , c — 1. 

-Tfc — Qi,i+fci fc > 1. 


We obtain: 


lim A • 1 = A 1 /, (3.3) 

Z—^OO 

y 0 = lim = T, (3.4) 

z —yoo 

Y\ = lim A~ 1 Qa + 1 = A” 1 Jr 0 + J, (3.5) 

2—^OO 

y fc = lim A^ 1 ^! i+fc-i = A _1 /r fe _i, k >2, (3.6) 

Z—^OO 

and 

OO 

We = e. (3.7) 

k=0 


It can be easily verified that the Markov chain > 0} satisfies Definition 1.1. Hence, it 

belongs to the class of asymptotically quasi-Toeplitz Markov chain. From (13.4D - (13.61) . we get 

Y{z) = T + Iz + &T x lT(z)z, \z\ < 1, (3.8) 


where 

TT))m' 
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here 


= 


O, 

Irwv © S®\ 

T 0 © D 0 © E 0 © S® l , 

Zg-l(z), 

®l(z), 


l' < l - 1,1 = 2,c, 
l' = l — 1, l = 1, c, 
l' = 1,1 = 0,5-1, 


l' = 9,1 = 0,5 - 1 , 

l' = 1,1 = g^c, 


Ir © (D r © £? r ) © I M i © <j® r , l' = l + r, r = 1 ,5 — 1 — l, l = 0 ,5 — 2, 
/rVF © E r © I M i © ? 0r , 


l' = l + r, r = g + 1 — l,c— 1 — 1,1 = 0 ,5 — 1 , 


iijw © E r © © ? 0r , 

ViW. 


I' = l + r, r = 1 , c — 1 — l, l = g, c — 2 , 


V = c, l = 0, c — 1, 


,,0) = Ir© 


71—1 


I W ®E n + Z n D(z) - ^ D k Z k © IyMa-n 


k =0 




©„(z) = T 0 © D( 2 ) © ^0 © S® n , n = g,c — 1, 


®c(z) =T 0 ® D(z) © E(z) © 5® c , 


*n(z) = Irw © z~ n (E(z ) - ]T E k z k ) © / M c-„ © ? 0n , 

V fc=0 / 


71 = 1, C. 


Substituting the expressions of the matrices T, /, A and T(z) in to (13.81) . we see that the matrix 
function r(z) is of the form 


Y(z) = 


Y u Y\ 2 
Y 22 (z) 


(3.9) 


where 


9-2 


( 


Y n = : 


9-3 

9-2 


Orwv Irwv © ? 


Irwvm © ? 


Irwvms-3 © ? 

O 


RWVMa- 2 J rwv ^-2 


M' 
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0 


/ 


Vl2 = 


9~2 \ IrwVMi-* ® ? ) v? -a 


RWVJ2IZ o M^x wv Ei= ff -1 Mi 


12 


M*) = 


OrwVMs-i Y‘22 

^22 ( z ) ^22 ( z ) 


V 12 - 
1 22 ~ 


l RWVM9 - 1 ® ^)O flW - V r Af9- 1 xi?wyx;= =(9+ i)M") > 


= (^A- 1 


® ^o* J i xiiivy E^(, + i) 



( A- 1 © g (z) + I 

A gg ^g,g+ 1 

Agg 1 r 9 ,c-l 

A gg ^c-g(z) 

II 

aa 

AJ+i.a+i ’ Irwv ® H® 9+1 

Ag+l,g+l®9+l( z ) + ^ 

a 9 -( 2 , 9+2 'W®< s+2 ■ 

Ag +1 ,g +i rg+l,c-l 

^g+2, 9 +2^9 + 2.c-l 

^g+l,g+l'F-9- 1 ( 2 ) 

V+2, 9 +2 l ^ c -9- 2 ( 2 ) 




A c -i,c-i® c_1 ( 2 ) + 1 

Ai-u-ifiW 


\ 


A c,l ■ Irwv ® 

Ac,c©c(z) + I 


\ 


/ 


here 


A/ ; =diag{(T r ,r = 1, i?} © w = 1,W} 

© diagjA^, v = 1, F} © [dia^{s m , m = 1, M}]® z , l = g, c, 


IV = Irvk E r © I M 1 © V, Z' = l + r, r = 1, c — 1 - l, l = g, c - 2. 


It is seen from (13.91) that the matrix Y(z) is reducible and the matrix F(l) contains only one 
irreducible stochastic diagonal block Y 22 (1)■ According to the Proposition 1.2, we get a sufficient 
condition for ergodicity of the Markov chain {£i,t > 0} which is the fulfillment of the following 
inequality: 


[det(zl - Y 2 2 {z))]' | z= i > 0. 


(3.10) 


A more constructive form of ergodicity condition (13.101) in our system is given by the following 
theorem. 


13 














Theorem 3.1. A sufficient condition for ergodicity of the Markov chain > 0} ; is the fulfill¬ 
ment of the inequality 


A + ^<i, 

hi h2 


where 


hi — e M9-!i h2 — ^2*->o e M c—1 5 


and X,;, i = 1,2, are the unique solutions of the following system of linear algebraic equations: 

I X 1 [5®» + S® s (/ ms -i®?)]=0, 

1 X ie = l, 


X 2 [5® c + Sq C (I M c~i <g) Q] = 0, 

X 2 e = 1. 


Proof. The proof of this theorem follows the steps given in the paper 
consequently we omit the details here. 


Breuer et al. 


( 20023 ) and 

□ 


Next, we assume that the stability condition is fulfilled. Denote the steady state probabilities 
of the Markov chain as 


P ^o, b, v, 7, , • • • , 

= lim P < ot = o, bt = b, v t = v, 7 1 = 7 , = m^ l \ • • • , m+ bt ' > = m^ 

t —>00 L 

Let the row-vector Pj of dimension K denote the stationary-state probabilities 


Pi = P (i, b, v, 7 , m^\ ■ ■ ■ , i> 0, 

and define the infinite-dimensional probability vector P = (Pi,P 2 ,---). When the system is 
stable, P is the unique solution to PQ = 0 and Pe = 1. However, it is still an open problem to 
express P in the closed form (e.g., generating function), and even though Q is highly structured, 
P cannot be expressed in a tractable analytical form. In the following, we will adopt the algorithm 
of asymptotically quasi-Toeplitz Markov chains to solve the equilibrium equation PQ = 0. 
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4. Algorithm 


In this section, we use the numerically stable algorithm which has been elaborated in Klimcnok & Dudin 


(2006) to calculate the stationary distribution of the system. The algorithm is based on censoring 


technique of the asymptotically quasi-Toeplitz Markov chain and consists of some essential steps 
which are given in the following Algorithm 4.1. 

Algorithm 4.1. Computing the probability vectors Pj.‘ 

Step 1: Calculate the matrix G as the minimal nonnegative solution to the non-linear matrix 
equation: 


G = Y(G). 


(4.1) 


We should note that the first RWV rows of the matrix G coincide with corresponding rows of 
the the matrix T(l) and the rest RWV V',, M n entr ies of the matrix G can be calculated by the 


iterative method which is available in 


L uc a nto n i (199%). 


Step 2: Chose the integer ko as a minimal value of k for which the norm ||G — Y^=k Qk+i,nG n ~ k \ 
of the residual is less than the preassigned value e. 

Step 3: Calculate the matrices Gk 0 - 1 , Gk 0 ~2 , ■ ■ ■, Go from the following recursive equation 


/ 00 \ 1 

Gk — J ^ ( Qk+l,n,Gj — lGj—2 ’ ’ ’ Gk+x I Qk+l,ki 
V j=k +1 / 

k = 0, 1, • ■ • , fco — 1, with the boundary condition Gk = G, k > ko. 

Step f: Calculate the matrices H^j, j > i, i > 0,using the formidae 

OO 

l It .j — H T *y ( Hi,kGk—iGk—2 ■ ■ ■ Gj. j (■'’ i . i > 0, 

k=j+l 

where Gk = G, k > ko. 

Step 5: Calculate the matrices Fj, j > 0, according to the recursive relation 


(4.2) 


(4.3) 


j -1 


Fn = 


i, /•) ^ n,j) j > 1 . 


(4.4) 


i=0 


Step 6: Calculate the vector Po as the unique solution to the system of the linear equations 

Po(-tfo,o) = 0, 

OO 

P 0 E Fj e = 1. 

j =0 
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Step 7: Calculate the vectors Pj, j >1, as 

Pj = P 0 F j: j > 1. 

Remark 4.1. We can eliminate the infinite sums in Step 6 as the following method: pick N > k$, 

N 

and replace the second equation in Step 6 with P 0 Fje = 1, such that the inequality P^e < eo 

3=0 

is satisfied, where eo is some preassigned small number. Note that the inverse matrices that appear 
in the algorithm exist and are nonnegative. So the computation of the algor ithm is st able. The 
algorithm can be realized in framework of extension of U SIRIUS++” software 


Dudin et al. 


200i) 


that was developed using the results of the asymptotically quasi- Toeplitz Markov chains. As such 
the determination of the stationary distribution is straightforward. 

5. performance measures 


Having the stationary distribution P,, i > 0 been calculated, we are able to calculate some 
stationary performance measures of the system under consideration in this section. 

Let P i — (P jo, P , , P if) , ^ 0, where P u, of size 1 x RWVM b is a row vector, 0 < b < c. 
Denote P(i, b) the joint probability to have i repeated customers in the orbit and b customers on 
the servers at arbitrary time. From the results of the above section, we get 

R W V M M 


P{h b) = P lb e = X! P 


r= 1 v=\ 7=1 m (l) = l 




(1) 


, m 


{b) 


= P*e 


1 - M b 1 - M b+1 

WV- - 7 ^ + 1 ,wv 


1-M 


1 - M 


where e 


RWV^. + 1,RWV^1 


of size K x 1 is a column vector of suitable size having 


ones as its RWV \jf~h + 1th to RWV 1 th entries and zeros as the rest entries. 

Corollary 5.1. 


1. The probability that there are i repeated customers in the orbit at arbitrary time 


p (l # ) = Y p ( i,b ^ = p * e > 1 - °- 


6=0 


2. The probability that there are b busy servers at arbitrary time 

OO 

P (*’ 6 ) = Y P ^' b ^ h = 27c- 


i =0 
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3. The mean number of busy servers 

C 

L b = Y^bP(*,b). 

6=1 

4■ The mean number of repeated customers in the orbit 

OO 

Lorb = ^ ^ 
i =1 

5. The mean number of customers (primary customers and priority customers) in the system 

L s — L or b + L\). 


6. The blocking probability for an arbitrary primary customer 


9 oo 


Pb i = l- — ~ n ">Tj( k ~ n )i.lR®D k ® I VMg -n)e. 


n= 1 i =0 


k =0 


7. The blocking probability for an arbitrary batch of primary customers 


1 


9 oo 


Pwi = i - t— TTTp ^' 9 ~ n ">Tj( lR ® Dk ® 


n= 1 i =0 k =1 

8. The blocking probability for an arbitrary priority customer 

1 


p b2 = 1 - P ^ i ' c ~ n ) - n )(^w ® E k <g> I M c—n )e. 


7i=l i=0 


k =0 


5. The blocking probability for an arbitrary batch of priority customer 

1 c oo n 

Pbb2 = 1 - t— Y ~]P(i,c — ri) y^flRw ®E k ® I M c-n)e. 

62 n = l i=0 k= 1 

1 0. Busy period: If we define the busy period B of this queueing system with repeated customers as 
the period that starts at the epoch when an arriving batch of customers (primary or priority) 
finds an empty system (all servers are idle and no repeated customer in the orbit) and ends 
at the departure epoch at which the system is empty again, then the mean busy period is given 
by 


E(B) = 


1 


Ai + A 2 \P( 0,0) 


1 


- 1 


Moreover, using the approach based on Markov renewal processes, we can obtain the following 
statement. 
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Proposition 5.1. The generating function 11 ( 2 ;) of the queue state stationary distribution at the 
departure epochs is characterized by the following functional-differential equation: 


n'(*) = n (z)A( z ), 


where 

A(z ) = B~\z)B'(z) 


B 1 (z)A(T 1 1 0 L B{z) + z i (T 1 1 <g> I wv j 2 c k=0 M k )B{z), 


- 1-1 


-1 


-l/rj-i-l 


B(z) = [T -zl- A _1 Jzr(z)][T - zl - A-ffzT(z)]- 1 . 


-1; 


1-1 


Proof. The proof can be implemented similarly as the proof of 
here. 


Breueretal 


( 2002 ) and is omitted 

□ 


6. Optimization 


In this section, we borrow some ideas and methods from 


Trivedi et al 


( 2000 ) to investigate 


the optimization problem that how to choose the optimal values of g and c so as to minimize 
the priority customer blocking probability as well as the primary customer blocking probability. 
Therefore, we obtain a multi-objective optimization problem in which the number of total servers c 
and the parameter g which is equal to the number of non-guard servers are the decision variables. 
However, due to the complexity, we may fix c and only consider g as the decision variable to simply 
the optimization problem. In the case of the two objectives, we can choose either P^i or ff /2 as 
the objective function to be minimized and impose a constraint on the other one. If we regard 
the two blocking probabilities as two binary functions of g and c two variables, then we write two 
representative optimization problems as follows: 


(I) Given T*., D^^E^, k = 0,1, • • • ,5 and c determine the optimal integer values of g such that 

Mininize Pb\ ( g) 
subject to Pb 2 (g)<Po- 

(II) Given T)., D^^E^, k = 0,1, • • • , S determine the optimal integer values of g and c such that 

Mininize c 

subject to Pbi{g,c)<pi, and Pb 2 (g, c) < p- 2 . 
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Here, the constants po, p\ and p 2 are three preassigned numbers. 

We first consider the optimization problem (I). It is clear that Piaig) > Pb 2 {<] — 1) (see Fig. 3). 
Thus, we can get the largest value of g such that P^ig) < po- Then based on the monotonicity 
property Pbi(g) > Pt,\ (.9 + 1), we know that such a value of g will minimize Pbi(g)- Therefore, we 
can obtain the optimal value of g using a simple one dimensional search over the range { 1 , • • • , c— 1 } 
for g such that g* = maxjglPft/g) < Po}- 

Then we study the optimization problem (II). Since we cannot derive the analytical solution 
of this system, in the next section, we will resort to the numerical examples to illustrate the 
optimization problem (II). Here, we provide an algorithm to obtain the minimum value of c. 

Algorithm 6.1. Solving the optimization problem (II) 

Step 1: Set the values p\ and p 2 and let c=2; 

Step 2: For g = 1 to c— 1, find the feasible region for g: 

© = {g\Pbi(g,c) <pi, and P b 2 {g,c) <p 2 }; 

Step 3: 7/0/0, then stop and obtain c* = c. Otherwise, put c = c + 1 and goto step 2; 

Step 4: Obtain c* = c. 

7. Application to cellular wireless networks and numerical examples 

In this section, we present an application to cellular wireless networks for the model under 
discussion and some numerical results that show the effect of the system parameters on the per¬ 
formance measures and channel allocation schemes in the cellular mobile networks. In the cellular 
mobile wireless network, the regional service area is divided into multiple adjacent cells, each of 
them served by a base station (BS) with a limited number of channels. Here, we assume the number 
of channels is equal to c, c > 2. A mobile subscriber who has wireless terminal such as mobile 
phone for voice service, wireless data terminal for data service, dual terminal for voice and data, 
and so on, communicates via radio links to a BS, one for each cell. Generally, there are two types of 
calls in a specific cell in the cellular mobile networks: originating calls, i.e., calls originating in that 
cell, and handoff calls, i.e., current ongoing calls caused by mobile users moving from one cell to 
that cell. We assume that the originating calls arrive at the system according to a BMAP with the 
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characterizing matrix sequence ( D n : n € No) and the handoff calls arrive at the system according 
to a BMAP with the characterizing matrix sequence (E n : n £ No). When a batch of new calls 
originates in that cell, some idle channels are required to the BS for the new calls setup. If there 
are idle channels in the BS of that cell, then the BS assigns some idle channels to the originating 
calls; otherwise some new calls (maybe all the batch) in that batch are blocked and the blocked 
mobile subscribers may try their luck later. Similarly, if some idle channels are available in the BS 
of that cell, the batch of handoff calls is successfully handed over without an interruption; otherwise 
some handoff calls (or all the batch) are blocked and retry for service after a short period. Any 
blocked call (originating or handoff) retries for service where the rates of the retrials depend on a 
changing environment. Thus, we use the MMPP with the characterizing matrices (To, Ti) to model 
the retrial process. The quality of service is disturbed seriously if a handoff call becomes blocked 


when it crosses the cell boundary. Generally, the base station gives priority to a handoff ca,. 
an originatin g call. The usual way to do this is the scheme with guard channels, see 


Guerin 


and 


Hong & Rappaport 


1 over 


19SS0 


(1996i). Here, we assume the number of non-guard channels is equal to g. 


Thus, the number of guard channels is equal to c — g. Furthermore, we suppose the conversation 
time of an ongoing call (originating or handoff) by a channel is of PH type with a representation 
(g,S). Then the cellular mobile wireless network can be modeled as a BMAPi,BMAP 2 /PH/g,c 
retrial queueing system where the rate of individual repeated attempts from the orbit is modulated 
according to a MMPP. The originating calls can be regarded as primary customers and the handoff 
calls can be regarded as priority customers. 


In the following, we give some simple numerical examples that illustrate the effects of various 
parameters on the system performance measures of the cellular mobile wireless network. The algo¬ 
rithms 4.1 and 6.1 developed in sections 4 and 6 have been written into some MATLAB programs. 
Realization of these algorithms on computer does not meet any difficulty. Only the calculation 
time can be long especially when the dimension of the system is large. For example, if we let c = 8, 
then the dimension K = 4088, and the calculation time is 98s. For c = 100, the dimension is 
equal to 1.0141e+031, and the calculation time is very long. Thus, the calculation time increases 
drastically with the increase of the number of servers c. Note that in all the below examples, we 
choose the parametric values in a way such that the system is stable. Numerical results are showed 
in Tabs. 1-3 and Figs. 1-3. For comparison, in each of the pictures, we also plot three curves which 
correspond to cr = 6.6429, 13.2857 and 26.5714. 
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For the convenience of the numerical calculation, here, we suppose that handoff calls and 
originating calls arrive to the system according to two independent special BMAPs (In fact, they 
are MAPs) which are characterized by the following matrices: 


Do — A o 


-11 

5 


2 


-20 


D\ — \ Q 


8 1 
3 12 


Dk = O 2 , k > 2 , 


Eq = A h 


-3 

1 


0 

-2 


E\ = A h 


1 2 
0 1 


Ek = O 2 , k > 2, 


This means the stationary arrival rates of the originating calls and the handoff calls are Ai = 
10.6364A o and A 2 = 1.6667A/j, respectively. 


The repeated calls in the orbit repeat there attempts to reach a server according to the MMPP 
which is defined by the matrices: 


Tn — \r 


-15 3 

4 -19 


T\ — A r 


12 0 
0 15 


Thus, the mean retrial rate for the mobile engaged in the cell is a = 13.2857A r 
The PH service process is defined by the matrices: 


5 = 


-23 9 


? = (0.4, 0.6). 


14 -17 

So that the mean rate of service is // = 8.1288. 

From the above descriptions, we have W = 2, V = 2, M = 2 and R = 2. 


The joint distribution of the number of busy channels and the number of repeated calls in the 
orbit is presented in Tab. 1 with (c, g, A 0 , A/J = (8,6,2,2). From Tab. 1, we observe that the joint 
probability that there are three busy channels and zero repeated call at arbitrary time is equal to 
0.2148 which is the maximum value and the probability that there is no repeated call in the orbit is 
equal to 0.9084 which is the maximum value of the marginal probabilities (lateral). The probability 
that there are three busy channels at arbitrary time is equal to 0.2206 which is the maximum value 
of the marginal probabilities (longitudinal). 
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Tab. 1 The stationary join distribution of the system with (c, g , A 0 , A/,., A r ) = (8, 6, 2, 2, 2) 


i\b 0 

1 

2 

3 

4 

5 

6 

7 

8 

sum 

0 

0.0467 

0.1413 

0.2136 

0.2148 

0.1604 

0.0923 

0.0376 

0.0016 

0.0001 

0.9084 

1 

0.0001 

0.0006 

0.0020 

0.0047 

0.0091 

0.0149 

0.0208 

0.0013 

0.0001 

0.0536 

2 

0.0000 

0.0000 

0.0002 

0.0008 

0.0023 

0.0054 

0.0109 

0.0008 

0.0000 

0.0206 

3 

0.0000 

0.0000 

0.0000 

0.0002 

0.0007 

0.0022 

0.0055 

0.0005 

0.0000 

0.0092 

4 

0.0000 

0.0000 

0.0000 

0.0001 

0.0002 

0.0009 

0.0028 

0.0003 

0.0000 

0.0043 

5 

0.0000 

0.0000 

0.0000 

0.0000 

0.0001 

0.0004 

0.0014 

0.0001 

0.0000 

0.0020 


1 . 0 e -003 x 










6 

0.0000 

0.0001 

0.0007 

0.0055 

0.0334 

0.1661 

0.6986 

0.0667 

0.0051 

0.9761 


1 . 0 e -003 x 










7 

0.0000 

0.0000 

0.0002 

0.0019 

0.0130 

0.0731 

0.3474 

0.0337 

0.0027 

0.4721 


1 . 0 e - 003 x 










8 

0.0000 

0.0000 

0.0001 

0.0007 

0.0052 

0.0325 

0.1724 

0.0169 

0.0014 

0.2291 


1 . 0 e - 004 x 










9 

0.0000 

0.0000 

0.0002 

0.0025 

0.0211 

0.1460 

0.8533 

0.0846 

0.0069 

1.1145 


1 . 0 e - 004 x 










10 

0.0000 

0.0000 

0.0001 

0.0009 

0.0087 

0.0660 

0.4217 

0.0421 

0.0035 

0.5429 

sum. 0.0468 

0.1419 

0.2158 

0.2206 

0.1728 

0.1162 

0.0800 

0.0047 

0.0002 

0.999 

Tab. 2 The the optimal value g* for the optimization problem (I) with ( c , X 0l po) = (20,10,0.0001) 


\ Xh 1 2 

3 

4 5 

6 7 

8 9 

10 

11 12 

13 14 

15 

20 

1 

18 18 

18 

17 17 

17 16 

16 16 

16 

15 15 

15 14 

14 

13 

10 

18 18 

18 

17 17 

17 16 

16 16 

16 

15 15 

15 14 

14 

13 

20 

18 18 

18 

17 17 

17 16 

16 16 

16 

15 15 

15 14 

14 

13 

Tab. 3 The optimal value ( 

’* for the optimization problem (II) with (X r ,Pi,P 2 ) = 

= ( 10 , 0 . 001 , 0 . 0001 ) 


\ A 0 1 2 

3 

4 5 

6 7 

8 9 

10 

11 12 

13 14 

15 

20 

1 

8 11 

14 

16 18 

20 22 

24 26 

28 

29 31 

33 35 

37 

45 

5 

10 13 

15 

17 19 

21 23 

25 27 

29 

30 32 

34 36 

38 

46 

10 

13 15 

17 

19 21 

23 25 

27 29 

31 

32 34 

36 37 

39 

47 
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Fig. 1 L or b as function of the parameter g with (c, A 0 , A h) = (10, 2, 2). 



Fig. 2 Dependence of the blocking probability for originating calls on the value g with (c, A 0 , A h) 
( 10 , 2 , 2 ). 
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Fig. 3 Dependence of the blocking probability for handoff calls on the value g with (c, X 0 ,Xh) = 
( 10 , 2 , 2 ). 

The effects of the retrial rate and the arrival rate of the originating calls on the optimal value 
of g for the optimization problem (I) with the set of parameters (c,Xh,0o) = (20,10,0.0001) are 
reported in Tab. 2. We observe that the optimal value g* decreases monotonously as the arrival 
rate of the handoff calls increases. From this, we can predict that the optimal value g* increases 
monotonously as the arrival rate of the originating calls increases. However, the retrial rate of the 
repeated calls in the orbit has no significant impact on the optimal value g*. The appearance of 
this phenomenon may be based on the fact that the mean number of repeated calls is very small 
compared to the mean number of arriving calls. Hence it may be concluded that the arrival rates 
of the two types of calls are the major factors in deciding the optimal numbers of guard channels. 

Tab. 3 shows the effects of the arrival rates of the originating calls and the handoff calls on the 
optimal value c* for the optimization problem (II), where we set (\ r ,Pi,P 2 ) = (10,0.001,0.0001). 
As is to be expected, when the arrival rates of the two types of calls increase, the optimal value c* 
increases. Moreover, we observe that the optimal value c* increases faster with the arrival rates of 
the originating calls than that with the arrival rates of the handoff calls. Therefore, we understand 
that the arrival of originating calls plays a most significant role in deciding the optimal value c*. 
In other word, this type of calls whose arrival rate is largest will most strongly affect the optimal 
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design of the system. 

Fig. 1 illustrates the dependence of the mean number of repeated calls in the orbit on the 
value of parameter g where we set (c, A 0 , A h) = (10, 2, 2). As it can be seen from Fig. 1 that L or f, 
decreases monotonously as the value g increases, especially when g < 6. But when g > 7, the 
curves of the values Lori, are almost overlapping with x-axis. That is to say, there are very few 
repeated customers in the orbit when the number of guard channels is small. 

Figs. 2 and 3 show the dependence of the blocking probability for an arbitrary originating 
call and handoff call on the value of parameter g respectively when (c, \ 0 ,\h) = (10,2,2). From 
Figs. 2 and 3, we find that the blocking probability for an arbitrary originating call Pb\ decreases 
monotonously, but the blocking probability for an arbitrary handoff call Pb 2 increases monotonously 
as the value g increases, which agrees with the intuitive expectations. In addition, when g < 6, the 
curves of the values are almost overlapping with x-axis. It means that the blocking probability 
for an arbitrary handoff call P \,2 is roughly equal to zero when the number of guard channels is 
large. 

We must note that, in each of the pictures, the three curves are almost overlapping. This 
phenomenon tells us that the retrial rate a has little impact on the system performance measures. 

From what has been described above, we can draw the conclusion that the performance measures 
are mainly affected by the two types of calls’ arrivals and service patterns, but the retrial rate plays 
a less crucial role. 


8. Conclusions 


In this paper, we have investigated the BMAP/PH /N type retrial que ue wi th two types of cus- 


tomers. Our work can be considered as an extension of 


Breuer et a! 


(2003) and 


Dudin &; Klimcnok 


(2012). The behavior of this model is described by a multi-dimensional continuous-time asymp¬ 


totically quasi-Toeplitz Markov chain. Sufficient condition for the ergodicity of the Markov chain 
is given. An algorithm for computing the stationary distribution is presented. Expressions for 
calculation of main performance measures of the model are derived. The optimization problem 
how to choose the optimal values of g and c is discussed and an algorithm is also provided. An 
application of the model to the cellular wireless network is implemented. The dependences of the 
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main performance characteristics of the model on its parameters are graphically demonstrated. We 
obtain the results that the retrial pattern has very little impact on the optimal values of guard 
channels and total channels, the number of busy channels and the blocking probabilities for the 
two types of calls. 

Future work will consider the BMAP/PH/c retrial queues with vacations or working vacations 
and their applications to computer communication networks. 
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